package edu.utah.seq.analysis;

import java.io.*;
import java.util.*;

import edu.utah.seq.data.*;
import trans.tpmap.WindowMaker;
import util.gen.*;

/**Use to find the peaks within EnrichedRegions*/
public class SubPeakPicker {
	//fields
	private EnrichedRegion[] enrichedRegions;
	private File[] treatmentPointDirs;
	private File[] controlPointDirs;
	private double scalarTreatment;
	private double scalarControl;
	private HashMap<String, PointData[]> treatments;
	private HashMap<String, PointData[]> controls;
	private int windowSize = 50;
	private int halfPeakShift;
	

	//constructors
	/**Picks the best peak based on rescanning each enriched region using the defined window size.*/
	public SubPeakPicker (EnrichedRegion[] enrichedRegions, File[] treatmentPointDirs, File[] controlPointDirs, int windowSize, int halfPeakShift){
		this.enrichedRegions = enrichedRegions;
		this.treatmentPointDirs = treatmentPointDirs;
		this.controlPointDirs = controlPointDirs;
		this.windowSize = windowSize;
		this.halfPeakShift = halfPeakShift;
		pickPeaks();
	}
	
	
	//methods
	private void fetchData(){
		System.out.println("\n\tLoading PointData...");
		treatments = PointData.fetchStrandedCombinePointData(treatmentPointDirs);
		double numberTreatmentObservations = PointData.totalObservationsMultiPointData(treatments);
		//System.out.println("\t\t"+(int)numberTreatmentObservations+" Treatment Observations");
		controls = PointData.fetchStrandedCombinePointData(controlPointDirs);
		double numberControlObservations = PointData.totalObservationsMultiPointData(controls);
		//System.out.println("\t\t"+(int)numberControlObservations+" Control Observations");
		//make normDiff scalars 
		scalarTreatment = new Double (numberControlObservations/numberTreatmentObservations).floatValue();
		scalarControl = new Double (numberTreatmentObservations/numberControlObservations).floatValue();
	}

	/**Main method of class. This method will then re window scan each EnrichedRegion and set the subWindows and bestSubWindow.*/
	public void pickPeaks(){

		//fetch data 
		fetchData();

		//window score each EnrichedRegion
		scoreTreatmentsControls();
		
		//find and assign best peak for each EnrichedRegion
		findPeaks();

	}

	/**Finds the best peak within each loaded enriched region.  Call scoreTreatmentsControls() or scoreTreatments() first.*/
	public void findPeaks(){
		for (int x=0; x< enrichedRegions.length; x++){
			SmoothingWindow bestSubWindow = findHighestPeak(enrichedRegions[x].getSubWindows(), 0);
			enrichedRegions[x].setBestSubWindow(bestSubWindow);
		}
	}
	
	/**Find the window with the largest score.  If multiple windows have the same best score then the start of the first and end
	 * of the last are used to create a new window and are returned.*/
	public static SmoothingWindow findHighestPeak(SmoothingWindow[] sm, int scoreIndex){
		if (sm == null || sm.length ==0) return null;
		//find highest score
		float bestScore = sm[0].getScores()[scoreIndex];
		for (int i=1; i< sm.length; i++){
			float testScore = sm[i].getScores()[scoreIndex];
			if (testScore > bestScore) bestScore = testScore;
		}
		//find smallest start and biggest end
		int start = sm[sm.length-1].getStart();
		int end = 0;
		for (int i=0; i< sm.length; i++){
			float testScore = sm[i].getScores()[scoreIndex];
			if (testScore == bestScore) {
				//set start?
				if (sm[i].getStart() < start) start = sm[i].getStart();
				//set end?
				if (sm[i].getStop() > end) end = sm[i].getStop();
			}
		}
		return new SmoothingWindow(start, end, new float[]{bestScore});
	}

	public void scoreTreatmentsControls(){
		//make window maker
		WindowMaker wm = new WindowMaker(windowSize,1);

		//fetch PointData
		PointData treatPD = null;
		PointData ctrlPD = null;
		int[] treatPos = null;
		int[] ctrlPos = null;
		
		//for each enriched region, these are sorted by chromosome
		String currentChromosome = "";
		for (int x=0; x< enrichedRegions.length; x++){
			//check chromosome and load PointData
			String testChromosome = enrichedRegions[x].getChromosome();
			if (currentChromosome.equals(testChromosome) == false){
				currentChromosome = testChromosome;
				//fetch PointData
				PointData[] pdT = treatments.get(currentChromosome);
				PointData[] pdC = controls.get(currentChromosome);
				//shift				
				pdT[0].shiftPositions(halfPeakShift);
				pdT[0].stripScores();
				pdT[1].shiftPositions(halfPeakShift);
				pdT[1].stripScores();
				pdC[0].shiftPositions(halfPeakShift);
				pdC[0].stripScores();
				pdC[1].shiftPositions(halfPeakShift);
				pdC[1].stripScores();
				//merge
				treatPD = PointData.combinePointData(pdT, true);
				ctrlPD = PointData.combinePointData(pdC, true);
				//get positions
				treatPos = treatPD.getPositions();
				ctrlPos = ctrlPD.getPositions();
			}
			//find indexes		
			int[] treatIndexes = treatPD.findIndexes(enrichedRegions[x].getStart(), enrichedRegions[x].getStop());
			int[] controlIndexes = ctrlPD.findIndexes(enrichedRegions[x].getStart(), enrichedRegions[x].getStop());

			//fetch positions
			ArrayList<Integer> positions = new ArrayList<Integer>();
			for (int i=treatIndexes[0]; i< treatIndexes[1]; i++) positions.add(new Integer(treatPos[i]));
			for (int i=controlIndexes[0]; i< controlIndexes[1]; i++) positions.add(new Integer(ctrlPos[i]));
			int[] pos = Misc.integerArrayListToIntArray(positions);
			if (pos == null){
				System.out.println(currentChromosome+"\t"+enrichedRegions[x].getStart()+"\t"+enrichedRegions[x].getStop());
				Misc.printExit("\nNo reads found within region? Are you sure you have entered the same PointData as what was used in ScanSeqs?\n");
			}
			Arrays.sort(pos);

			//create windows
			int[][] windows = wm.makeWindows(pos);
			//assign bp positions
			for (int i=0; i< windows.length; i++){
				windows[i][0] = pos[windows[i][0]];
				windows[i][1] = pos[windows[i][1]]+1;	//last base isn't included
			}

			//for each window
			ArrayList<SmoothingWindow> smAL = new ArrayList<SmoothingWindow>();
			for (int i=0; i< windows.length; i++){

				//fetch scores
				float tSum = treatPD.sumPositionBP(windows[i][0], windows[i][1]);
				//don't proceed if tSum ==0
				if (tSum == 0) continue;
				float cSum = ctrlPD.sumPositionBP(windows[i][0], windows[i][1]);

				//calculate score, save window
				float score = normalizedDifference(tSum,cSum);
				smAL.add(new SmoothingWindow (windows[i][0], windows[i][1], new float[]{score}));
			}
			//convert and add to enriched region
			SmoothingWindow[] sm = new SmoothingWindow[smAL.size()];
			smAL.toArray(sm);
			enrichedRegions[x].setSubWindows(sm);
		}

	}
	
	/**Calculates a normalized difference after scaling. Assumes tSum or cSum !=0*/
	private float normalizedDifference(float tSum, float cSum){
		double t;
		double c;
		//scale cSum?
		if (tSum == 0){
			//scale cSum
			t = tSum;
			c = cSum * scalarControl;
		}
		//scale tSum
		else {
			t= tSum * scalarTreatment;
			c = cSum;
		}
		return (float)((t-c)/Math.sqrt(t+c));
	}
	

	
}
